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Abstract. 

In a smooth flow, the leading-order response of trajectories to infinitesimal 
perturbations in their initial conditions is described by the finite-time Lyapunov 
exponents and associated characteristic directions of stretching. We give a description 
of the second-order response to perturbations in terms of Lagrangian derivatives of the 
exponents and characteristic directions. These derivatives are related to generalised 
Lyapunov exponents, which describe deformations of phase-space elements beyond 
ellipsoidal. When the flow is chaotic, care must be taken in evaluating the derivatives 
because of the exponential discrepancy in scale along the different characteristic 
directions. Two matrix decomposition methods are used to isolate the directions of 
stretching, the first appropriate in finding the asymptotic behaviour of the derivatives 
analytically, the second better suited to numerical evaluation. The derivatives are 
shown to satisfy differential constraints that are realised with exponential accuracy in 
time. With a suitable reinterpretation, the results of the paper are shown to apply to 
the Eulerian framework as well. 
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1. Introduction 



We consider a collection of coupled ordinary differential equations (ODEs) associated 
with a vector field [f] (a continuous-time dynamical system). The general solution to 
these equations defines a flow — a mapping from the phase-space domain onto itself .\ For 
smooth vector fields, the flow can be regarded as a smooth coordinate transformation 
(a diffeomorphism) from the set of initial conditions to the state of the system at some 
later time [2, p. 276]. The coordinates describing the initial conditions are called the 

\ In the mathematics literature the term flow is reserved for autonomous systems, the term 
transformation being preferred for the solutions generated by time-dependent vector fields [2, p. 96]. 
Because of the compelling nature of the fluid analogy, flow will be used in the more general sense in 
this paper. 
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Lagrangian coordinates, and those describing the state at a later time are called the 
Eulerian coordinates. 

Generally, even a smooth vector field will lead to chaotic dynamics, § with 
trajectories of nearby phase-space elements diverging rapidly from each other, at an 
average rate close to exponential for long times. The allowable rates are called the 
Lyapunov exponent of the flow [4], and are associated with characteristic directions 
of stretching. For chaotic flows, the transformation to Lagrangian coordinates becomes 
exceedingly contorted, and in practise it can no longer be inverted, due the exponentially 
growing errors on the position of phase-space particles. Nevertheless, the presence of 
chaos can actually be advantageous because it leads to a large separation of timescales 
along the different characteristic directions of the flow. This separation of scale was 
used by Boozer [5], Tang and Boozer [6,7], and Giona and Adrover [8,9] to study fluid 
mixing and the dynamo problem. 

Many equations of fluid dynamics are "advective" in nature, in that they describe 
the motion of a scalar or vector field as it is dragged by a flow, and possibly influenced by 
other effects such as diffusion and sources. Examples are the scalar advection-diffusion 
equation [10] and the induction equation for a magnetic field [11]. When expressed in 
Lagrangian coordinates, the advective term drops out of these types of equations. For 
scalar and vector advection-diffusion equations, one is left with a diffusion equation 
with anisotropic diffusivity. The anisotropic diffusivity arises because the Jacobian of 
the transformation between Lagrangian and Eulerian coordinates is not orthogonal. In 
a chaotic flow, the complexity of the transformation leads to a singular Jacobian, which 
is reflected in an exponentially-growing diffusivity that enhances mixing [6,12]. 

In Refs. [6, 12], the singular Jacobian is decomposed to isolate the dominant 
direction of enhanced diffusion, leading to an expression in terms of the finite-time 
Lyapunov exponents and characteristic directions of stretching. To get a full solution of 
the advection-diffusion problem, the Lagrangian derivatives of the finite-time Lyapunov 
exponents and of the characteristic directions are needed. This is because the Lagrangian 
coordinate frame is position-dependent, so when derivatives of vector fields are taken 
one must also differentiate the basis vectors themselves (this is the same procedure 
as in covariant differentiation, or when fictitious forces appear after transforming to a 
rotating frame). Thus the necessity of obtaining Lagrangian derivatives of the vectors 
defining the coordinate frame. Because Lagrangian coordinates are also stretched with 
respect to Eulerian coordinates, the derivatives of the characteristic rates of separations 
(as characterised by the finite-time Lyapunov exponents) are also needed. 

The problem of finding the asymptotic form of Lagrangian derivatives of the 
coordinate transformation induced by a flow has been addressed previously by Dressier 
and Farmer [13] and Taylor [14] in a different context. They examined the asymptotic 
behaviour of the Hessian, the quadratic form consisting of the second derivatives of the 
flow. The Hessian is the term that follows the Jacobian matrix in a Taylor expansion 

§ Important exceptions are ID vector fields, and 2D autonomous vector fields [3]. There are also other, 
more restricted classes of nonchaotic vector fields. 
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of the coordinate transformation form Lagrangian to Eulerian coordinates. Dressier 
and Farmer call the growth rates of the Hessian generalised or higher- order Lyapunov 
exponents. Their motivation lay in characterising the growth of nonlinear distortions 
of geometric quantities evolving under the influence of one-dimensional maps. The 
Lyapunov exponents quantify the leading order stretching of an infinitesimal ellipse 
moving with the phase fluid, and the generalised exponents describe deviations from 
an elliptical shape. Dressier and Farmer provided numerical results only for the largest 
Lyapunov exponent, because no numerical method exist to evolve the Hessian in a 
numerically stable manner that is not susceptible to limited precision. We provide such 
a method here. 

In the Eulerian picture, the higher-order exponents also characterise the growth 
of extrinsic curvature of curves and surfaces embedded in a flow (material lines and 
surfaces), and so can be connected to work describing the evolution of the curvature 
of material lines in turbulent flows that originated with Batchelor [15-19], and more 
recently was applied to chaotic flows [20-22]. Our results — derived in the Lagrangian 
frame — can easily be adapted to the Eulerian picture, and provide a comprehensive 
view of second-order deformation processes in chaotic flows. The connexion between 
the Eulerian and Lagrangian pictures is made in an appendix. 

Our approach is similar to Refs. [13, 14] but is more general and applies to flows 
rather than maps: we aim to give estimates of the asymptotic growth rates of the 
Lagrangian derivatives of finite-time Lyapunov exponents and characteristic directions 
by appealing to arguments of "genericity" of the quantities involved, thus showing that 
the estimates will hold in essiantially all cases. These arguments are formal in nature, 
in the sense that they do not provide a mathematical proof of the results, since there 
will always exist flows that can be specifically chosen to violate any of the estimates. 
In particular, there could be flows with degenerate Lyapunov exponents (finite-time 
Lyapunov exponents that are degenerate for short periods of time do not concern us). 
We expect that even for degenerate exponents the results of the paper are applicable 
in a limited form. In practise the asymptotic behaviour holds to great accuracy for all 
flows examined. This will be verified in the numerical section of the paper. 

To obtain the estimates of asymptotic behaviour of derivatives, we perform a 
singular value decomposition (SVD) of the tangent mapping of the flow, and differentiate 
the ODEs derived by Greene and Kim [23] directly. A careful analysis of the equations, 
with the assumption of a nondegenerate spectrum and a bounded attractor as in 
Goldhirsch et al [24], leads to our asymptotic forms. We find that the Lagrangian 
derivatives along diverging directions of the finite-time Lyapunov exponents grow 
exponentially at the characteristic rate of that direction. This is consistent with 
the intuitive notion that small displacements in those directions will be exponentially 
amplified, so one expects derivatives to grow as well. For contracting directions, the 
opposite is true: the Lagrangian derivatives converge exponentially to time-asymptotic 
values, but often do so at a slower rate than the characteristic rate. The Lagrangian 
derivatives of the characteristic directions of stretching have a more complicated 



Derivatives and Constraints in Chaotic Flows 



4 



behaviour, and it is not always true that a derivative along an expanding direction 
will diverge — the intuitive picture fails. 

The asymptotic behaviour derived using the SVD method can be used to recover so- 
called differential constraints on the finite-time Lyapunov exponents and characteristic 
directions. Such constraints were first derived in two dimensions by Tang and Boozer [6] 
and Giona and Adrover [8] and were later extended to three dimensions by Thiffeault 
and Boozer [25] . These earlier derivations provided limited results and were difficult to 
generalise to higher dimensions. In particular, the convergence rate of the constraints 
was not obtained, making it difficult to gauge their effect in approximations. The 
derivation given here overcomes these difficulties. 

Whilst the SVD method is transparent and useful for theoretical derivations and 
interpretation, it is not well-suited for numerical purposes [26]: it possesses troubling 
singularities and involves a needlessly large number of equations to evolve. A better 
method is the QR decomposition [24,26], also known as the continuous Gram-Schmidt 
orthonormalisation method because it is a time-continuous version of earlier methods 
that involved re-orthonormalising a set of vectors evolved using the tangent map of 
the flow [3,27,28]. As for the SVD method, we adapt the QR method to finding the 
Lagrangian derivatives of the finite-time Lyapunov exponents and of the characteristic 
eigenvectors. The QR method can be used to verify the differential constraints 
mentioned above, since they depend on delicate cancellations that require the high 
accuracy afforded by the method herein in order to be convincingly established. 

The outline of the paper is as follows. In section 2 we introduce the basic framework 
and notation necessary to the subsequent development. The central object of study is 
the metric tensor transformed to Lagrangian coordinates, and its diagonal form that 
contains all the information on the characteristic separations and directions. Section 3 
describes the direct method of evolving the Hessian, where its governing equations are 
obtained by a variation of the ODEs and integrated directly. This method is not very 
useful numerically because the exponential blowup of the elements of the Hessian leads 
to issues of limited precision, but illustrates the basic principles and can be used to 
check the results of more complex methods, for short times. 

A more powerful method is introduced in section 4, the SVD decomposition method. 
This allows derivation of the asymptotic behaviour of the Lagrangian derivatives. The 
asymptotic behaviour is exploited in section 5 to investigate the properties of the 
Hessian, and to give a new and powerful derivation of differential constraints in chaotic 
flows. 

In section 6 the QR decomposition is used to develop a suitable numerical method 
for evolving the various Lagrangian derivatives. The numerical method is then used 
to verify to high precision the geometrical constraints derived in section 5. Section 7 
consists of a brief summary of the main results of this paper and of possible future work, 
as well as a discussion of possible applications. 
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2. Characteristic Directions of Trajectory Separation 

We begin with a brief overview of the concepts and notation we shall use. We consider 
the n-dimensional dynamical system 

x = v(x,t) , (2.1) 

where the overdot indicates a time derivative, and v is a smooth function of x and t. 
The solution to (2.1) is a function x(t), with initial condition a; (to) — a - We can thus 
regard x(t) as a coordinate transformation from the set of initial conditions a to the 
state at time t; we write this transformation explicitly as x(a,t). || Following standard 
terminology, we call x the Eulerian coordinates and a the Lagrangian coordinates. 
The time-evolution of the Jacobian matrix M l p := dx % jdaf is given by 

n 

M* p = J2G\M e p , (2.2) 

where G % £ :— dv l /dx e ; the initial conditions are M % p = 5 1 p , since the coordinates x 
and a initially coincide. The nonsingular Jacobian matrix tells us how to transform 
vectors in a coordinates to vectors in x coordinates (M l p is the tangent mapping to the 
transformation x(a,t)). Now construct the matrix 

n 

g pq :=Y.M i p M l <l , (2.3) 
e=i 

called the metric tensor in Lagrangian coordinates or Cauchy-Green strain tensor [29]; it 
is symmetric and positive-definite, so it can be diagonalized with real positive eigenvalues 
and orthogonal eigenvectors and rewritten as 

n 

9pi = E A l OMp (2.4) 

(7=1 

with e CT (a, t) and A^(a, t) respectively the ath eigenvector of g pq (a, t) and corresponding 
eigenvalue. We refer to the A CT as the coefficients of expansion, because they represent 
the relative deformation of the principal axes of an infinitesimal ellipsoid carried by 
the flow. The coefficients of expansion can be used to define the finite-time Lyapunov 
exponents, 

X a (a,t) := l - log Mo, t). (2.5) 

The finite-time Lyapunov exponents, X a (a,t) describe the instantaneous average rate 
of exponential separation of neighbouring trajectories.^ The multiplicative ergodic 
theorem of Oseledec [30] implies that the infinite-time limit of A cr (a, t) exists and is 
independent of the initial condition a in a given ergodic domain, for almost all initial 
conditions. The infinite-time limit e^°(a) of the characteristic eigenvectors e a (a,t) also 
exists but depends on the initial condition. The Lyapunov exponents converge very 
slowly, whereas for nondegenerate exponents the characteristic eigenvectors converge 

|| To lighten the notation, we omit the explicit dependence of x on the initial time to- 
% We are implicitly assuming the Euclidean norm for vectors in Eulerian space. 
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exponentially fast [24] . The slow convergence of the Lyapunov exponents indicates that 
the instantaneous separation rate of neighbouring trajectories is not at all exponential 
on a typical attractor; only in the infinite-time limit do the trajectories show a mean 
exponential rate of separation. The instantaneous deviations from this exponential 
rate are very large. However, even though it may not be growing exponentially, an 
eigenvalue A CT associated with a positive Lyapunov exponent becomes very large after a 
relatively short time, and conversely an eigenvalue associated with a negative Lyapunov 
exponent becomes very small. It is thus an abuse of language, but a convenient one we 
shall use, to refer to the A CT 's as growing or shrinking exponentially. 

In this paper we shall assume that the eigenvalues A CT are nondegenerate and ordered 
such that A cr _ 1 > A CT . After allowing some time for chaotic behaviour to set in, we have 
that A (j A K for a < k. We shall make use of this ordering often in the subsequent 
development. 



3. Lagrangian Derivatives: Direct Method 

When the vector field v(x,t) is a known analytic function, explicit evolution equations 
for the spatial derivatives of the A^ and e^ can be derived, as pointed out in Refs. [6] 
and [31]. The method involves expressing the derivatives of A^ and e M in terms of 
derivatives of the metric tensor g. This is done by taking the Lagrangian derivative of 
the diagonal form (2.4) of g and dotting the resulting expression with the eigenvectors e M . 
We obtain 

^^Z^W).^- (») 

I? = (^) 

The derivatives of the metric are obtained from the Hessian K k r := d 2 x k /da q da r of the 
coordinate transformation via the relation 



dg 



obtained by differentiating the non-diagonal form (2.3) of g. The Hessian is symmetric 
in its lower indices, and it is computed by solving the evolution equation 

n n 

K k qr = £ G\ K\ T + £ X% M\ M' r , (3.3) 
t=i i,j=i 

where := d 2 v h /dx l dx j . For the existence and uniqueness of solutions to (3.3), 
it is necessary that d 2 v k /dx l dx j be Lipschitz, but it is sufficient that v be at least 
thrice-differentiable. Equation (3.3) is obtained by differentiating (2.2), and the initial 
condition is K = 0. The time derivatives (the overdots) are taken at constant a, so they 
can be commuted with Lagrangian derivatives. 

The linear part of (3.3) is the same as for (2.2), but now there is a nonlinear coupling 
term to M. Since for a chaotic flow the matrix M has at least one exponentially growing 
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eigenvalue, the elements of the Hessian K ctypically grow faster than M, owing to the 
nonlinear coupling. Numerically, the system (3.3) is thus unstable in an essential way, 
and we can expect the integration to give meaningless results (except along the dominant 
stretching direction) after the elements of K become too large. Nevertheless, if one is 
not interested in long-time behaviour the direct method can yield satisfactory results. It 
is not suitable for a detailed, accurate, long-time solution of the Lagrangian derivatives. 
We develop such methods in sections 4 and 6. 



4. Asymptotic Behaviour using the SVD Method 

4-1. Basic Method 

Any matrix, and in particular the Jacobian matrix M, can be decomposed into the 
product 

M = UFV T , (4.1) 

where U and V are orthogonal matrices and F is diagonal. The superscript T denotes a 
matrix transpose. This decomposition is called the singular value decomposition (SVD), 
and is unique up to permutations of rows and columns. The diagonal elements A a of F 
are called the singular values. Requiring that the singular values be ordered decreasing in 
size makes the decomposition unique (for nondegenerate eigenvalues). As can be seen by 
substitution of (4.1) into (2.3), the columns of V are eigenvectors of g, V qa = (e a ) q , with 
eigenvalues given by the diagonal elements of F T F, (F T F) aa = A 2 . The advantage 
of the SVD is that it separates neatly the parts of M that are growing or shrinking 
exponentially in size (as determined by the coefficients of expansion A a ). The SVD 
has the following interpretation: if we consider an infinitesimally small "ball" of initial 
condition obeying (2.1), it will deform into an ellipsoid under the action of the flow. 
The A a give the relative stretching of each principal axis of the ellipsoid, the orthogonal 
matrix V gives the principal axes of stretching in Lagrangian coordinate space, and 
the orthogonal matrix U gives the absolute orientation of the ellipse in Eulerian space. 
Constructing the metric tensor as in (2.3) eliminates the Eulerian orientation, retaining 
the essential features of the stretching. 

Greene and Kim [23] derived the equations satisfied by U, V, and F: 

A M = G^A„ (4.2) 



G^A 2 + Gy^A 2 ^ 



for /i^i/; 



for pi — v\ 



(v T v) 



A M A v 
'A 2 - A?/*"" 



Auv for p^; 



i 1 



(4.4) 
for n = v\ 



where 



G := U T GU, A:=G + G 7 
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Numerically, the SVD method has limitations: the number of quantities to evolve is 
large, and when A CT ~ A CT+1 the denominators become singular (See Refs. [23, 26] for a 
discussion of these issues). However, conceptually the method is very straightforward 
and transparent, and it gives explicit equations for all the quantities we are interested 
in, as opposed to the QR method (section 6) which needs to be corrected to yield the 
true value of the A CT . 

For large time, when A a /A K <C 1, o > k, we can approximate (4.3) and (4.4) by 

H<v; 

fi>u; 

\i = v. 



(U T U) 



G 


n < 




+G /Xl/ 


n > 







fi = 





(v T v) 



(4.5) 







where we have used 



' K/K 
A M /A, 

max 



A„ A 



n+i 



A, 



A, 



fi<v; 
\i = v. 



(4.6) 



The matrix 7 is symmetric, and defined such that (for nondegenerate A) all of its 
elements are decaying exponentially. Because of the ordering of the A CT , the diagonal 
element 7 W represents the element on the //th row or column that decays the slowest, 
that is, 7 W = max CT ^7 CT/1 . 

If we assume that the evolution of the system takes place on some bounded domain, 
then we know that G, and hence G and A, remains bounded (we also assume that the 
vector field v is smooth). Thus, the right-hand side of V T V in (4.5) goes to zero 
exponentially, and we can solve the equation perturbatively for large time, 



= v% + E v% / df + 0(7^) 

v J to 



for some constant V™, which can only be determined by solving the unapproximated 
equation (4.4). We conclude that for large t the matrix V has the form [24] 

V q , = ^V q , + V™, (4.7) 

Since V qjX = (e M ) 9 and 7 MM — > 0, the characteristic eigenvectors converge exponentially 
to their time-asymptotic value, (e^°) g = V™. The elements of the matrix U do not in 
general converge. The Lyapunov exponents have a very slow (logarithmic) convergence 
which does not concern us here, as we are considering timescales of fast (roughly 
exponential) convergence. 



4-2. Lagrangian Derivatives 

Having derived the equations of motion for the SVD of M, we can now take the 
Lagrangian derivative of these equations of motion, in a manner analogous to section 3. 
We define the quantities 



d log A u 



KflU 



E^ 



q,i 



qK * da* ' 



e 



KflU 



E^ 



v dVpa 
qK Vii da* ' 



(4.8) 
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which are simply the Lagrangian derivatives of A u , U iv , and V pv expressed in a convenient 
frame. Note that $ KMJ , = —^ KVIl and Q K/Jil/ = — Q KU ^- The Lagrangian derivatives 
may be regarded as the components fi of the curvature of the vector fields defined by the 
columns of V pu . From (4.2)-(4.4), we can find equations for the evolution of ^ KU , & KfMU , 



and Q Kfll/ , 



* w = E(^U *™ + E 4- + A « ( 4 - 9 ) 

cr cr 

^ = E [(^U ^ + (t/ T t>W - (t/ T t>W <iw 

cr 

+ E^^( C/T ^W (4-10) 
= E [(^U ev + (v T v) a , Q Kau - (y T v) av e WM 

cr 

+ E^^(^U (4-11) 



where 



X^:=E^^^^ (4-12) 

is symmetric in k and //. 

In Appendix A, the asymptotic behaviour of the derivatives is obtained from the 
equations of motion (4.9)-(4.11). We now summarise the main results of that section: 
for t ^> 1, by which we mean that the dynamical system has evolved long enough for 
the quantities A M to have reached a regime of quasi-exponential behaviour, we have that 
the Lagrangian derivatives defined by (4.8) evolve asymptotically as 

= max (A K , 7^) (4.13) 
= max (A K , 7/tK , 7w ) + ~max(A K ,l), (4.14) 

= max (7^A«, 7kk , 7wt , 7 W ) + e~„ ~ max ( 7/1I/ A« , 1) . (4.15) 

Recall that in all these cases the first index, k, denotes the characteristic 
eigendirection e K along which the Lagrangian derivative is evaluated. For a given k 
with A K 1, corresponding to an expanding direction of the flow, both $ KMJ , and \& KV 
grow exponentially with time (the constant is then irrelevant). Thus, Lagrangian 
derivatives of logA„ along an expanding direction k become more singular with time, 
to a degree commensurate with the separation of neighbouring initial conditions along 
that direction, as given by A K . 

Conversely, for k with A K < 1, corresponding to a contracting direction of the 
flow, both $ KMJ y and ^/ KU decrease exponentially with time, $ Ktu , converging toward zero 
and ty KU converging to a constant, . The convergence rate of both these quantities, 
however, is not necessarily equal to A K but may be slower (though still converging), as 
denoted by the max in (4.13) and (4.14). This slower convergence rate of \I/ has two 
sources: (i) From (4.7), the V K in the definition (4.8) of ^ KV converges at a rate 7kk ; (ii) 
The $ term in (4.9) limits the convergence to 7l/ „. 
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The interpretation of the long-time behaviour of Lagrangian derivatives of V, 
the characteristic eigenvectors, is less generic. For a contracting direction k, the 
derivatives @ Ktlv converge to the constant value at a rate of max(7 KK , 7^, 7^). This 
convergence rate is dominated by the convergence of the individual Vs in (4.8), as given 
by (4.7). For an expanding direction k, the specific behaviour of Q KjJlV depends on the 
relative magnitudes of the coefficients of expansion A. However, for a non-contracting 
direction k it is always true that 

(e KIW - e~j « A K , a k > 1, 

so the gradients of e along expanding directions grow much more slowly that those 
of log A. 

We close this section by considering the asymptotic behaviour of the Lagrangian 
derivative of the determinant, \g\, of the metric tensor g. From (2.4), 

\g\ = J2 A li 

V 

so that the Lagrangian derivative of the determinant along V K is 

^ ■= E V,* ^ log \g\ 1/2 = £ . (4.16) 

q v 

An equation of motion for Q K is obtained by summing (4.9) over u, yielding 

^ = E(^U V a + A K E X UUK , (4.17) 

a v 

where the $ term has dropped out. The term YLv X VVK is proportional to V • v, so for an 
incompressible flow the solution to (4.17) is Q K = 0. But in general, for a compressible 
flow with non-uniform V • v, we find 

fi K = max(A K , 7KK )fi K + l^°, t > 1, (4.18) 

where the time dependence of Q K is non-exponential. In contrast to the limiting 
rate 7 KK in (4.18) is due entirely to the convergence rate of V K) as given by (4.7). 



5. Properties of the Hessian and Constraints 

5.1. Symmetry of the Hessian 

We now make contact with the work of Dressier and Farmer [13] and Taylor [14] on the 
form of the generalised Lyapunov exponents, which describe the asymptotic behaviour 
of the Hessian. The Hessian, defined in section 3, can be recovered from the Lagrangian 
derivatives of section 4.2. Following Dressier and Farmer [13], we project the components 
of the Hessian onto the U and V bases and define 

K;„:=Y,U^Kl q V p ,V qu . (5.1) 

Writing K l vq = dM e p /da g , and using the SVD decomposition (4.1) for M, we find 

K K = A K ^ K 8 VK + A K + A u $^ K , . (5.2) 
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The asymptotic behaviour of the Hessian is easily derived from (4.13), (4.14), and (4.15), 

K^ = max(A K , k»K)K; v , (5.3) 

where K* v is a non-exponential function, as found for maps in Refs. [13,14]. 

Since the Hessian is symmetric in its lower indices, we could have equally well 
written 

K«„ = A K „„ 5^ + A K 0^ K + A M $„ KAt , (5.4) 

where we simply interchanged \i and v in (5.2). Equating (5.2) and (5.4), we find the 
relations 

A M (9^ + „„) = A„ fi^v, (5.5) 

A K {Q ilUK -Q u ^ K ) = A^$^ K - A M $^ K , /I, v, k differ. (5.6) 

Equation (5.5) defines n(n — 1) independent relations, whilst (5.6) defines n(n — l)(n — 
2)/2 relations. Thus, a total of n 2 (n — l)/2 quantities are dependent on the others and 
can be eliminated, which is exactly the number of dependent components of K*. The 
relation (5.6) can be solved for Q Kflv to yield 
1 



®KflV 2 




A„ aJ^ k+ a„ aJ^ 



(5.7) 



when //, u, and « differ. Similarly, (5.5) can trivially be solved for 9^. Equations (5.5) 
and (5.7) express all of the the G in terms of the \& and $. However, the 
asymptotic behaviour of Q KIIV cannot be recovered from these equations by simply 
substituting (4.13) and (4.14). The reason is that the asymptotic form (4.15) hinges on 
delicate cancellations between the $ and \l/ that are not manifest from simply looking 
at their equations of motion. For instance, in (5.7) the coefficient of each $ term grows 
exponentially, even though some of the 0's have been shown to converge. 

Although we are not using the SVD method to obtain numerical results, the 
considerations of this section also apply to the QR method of section 6. The 
relations (5.5) and (5.6) can then be used as a diagnostic tool to monitor the numerical 
results. 

5.2. Differential Constraints 

Rather than solving for the 0, if the flow is chaotic the relations (5.5) and (5.6) can be 
put to good use in another manner. The Lagrangian derivatives of U, as contained in $, 
are not quantities of great interest to us. They describe the sensitive dependence on 
initial conditions of the absolute orientation of phase-fluid elements in Eulerian space. 
This information is not necessary for solving problems in Lagrangian coordinates, and 
is too sensitive to initial conditions to be of use anyhow. We thus substitute the time- 
asymptotic form of $, given by (4.13), in the right-hand side of (5.5) and (5.6), yielding 

Qwii/ + = max , ^ 7^ , (5.8) 
©AH/* - ®vhk = t~~ max (A M , 7^) $ M!/K - ^ max (A,, , 7 MK ) , (5.9) 
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Table 1. The total number of type I constraints for low-dimensional systems, as given 
by (5.12). The rows denote n, the columns the number of contracting directions m < n. 
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where, as before, /x, z/, and k differ. Below we show that under certain conditions the 
right-hand side of (5.8) or (5.9) converges toward zero, giving us asymptotic differential 
constraints on \1> and 9. 

5.2.1. Type I Constraints For /i < z/, (5.8) can be written 

9^ + V V p = max (A„ , 7^) , n < v, (5.10) 

If the index v corresponds to a contracting direction (A„ 1), the right-hand side 
of (5.10) goes to zero exponentially fast, at a rate A u or 7^, whichever is slowest. In 
that case (5.10) is a constraint implying that for large t we have 

(6^ + 0, A,<A M , A,<1. (5.11) 

We refer to (5.11) as type I constraints. The total number of such constraints is 

Ni = m[n-±(m+ij\, (5.12) 

where n is the dimension of the space and m is the number of contracting directions 
(i.e., the number of negative Lyapunov exponents) possessed by the flow in a particular 
ergodic domain. Table 1 gives the number of type I constraints, N h as a function of n 
and m. 

In two dimensions, we typically have one contracting direction, so there is a single 
type I constraint. This is the same constraint that was derived in Refs. [6,25]. 

In three dimensions, for an autonomous flow, we typically also have one contracting 
direction. There are then two type I constraint. These constraints correspond to those 
derived in Ref. [25]. 

A special case of the type I constraints is obtained by setting v = n in (5.10), and 
then summing over /j, < n, to yield 

E 1^ (M 1/2 (en),) - E(en)^ log An ~ max (A n , ^) - 0. (5.13) 
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This constraint was discovered numerically and used by Tang and Boozer [6,7,31,32] 
and derived in three dimensions by Thiffeault and Boozer [25]. It has been used to 
study the anticorrelation between curvature and stretching [6,25,31] and to transform 
the advection-diffusion equation in Lagrangian coordinates to an approximate one- 
dimensional form [12]. (See also Appendix C for an Eulerian version of the constraint.) 
The present method not only gives the constraint in a very direct manner for any 
dimension n, but it also provides us with its asymptotic convergence rate, as determined 
by the right-hand side of (5.13). It also shows that the constraint (5.13) does not stand 
alone, but is the sum of several independent constraints. More details on the differences 
between this paper and earlier approaches are given in section 5.3. 



5.2.2. Type II Constraints Equation (5.9) implies that 

n n (KK K A M \ 

0,h/k - ~ max ^— — , — 7„ K , — 7 MK J . (5.14) 

We are interested in finding constraints analogous to the type I constraints 
of section 5.2.1. It is clear that unless both // and v are greater than k, the right- 
hand side of (5.14) is of order unity or greater, and so does not go to zero. We can 
assume without loss of generality that /i < v, so that 

Q^uk - ~ lp,K max (A„ , 7 MK ) , k < fx < u, (5.15) 

where we have used r ) VK <C 7 MK . Whether or not (5.15) is a constraint depends on the 
specific behaviour of A M A^/A K . Clearly, we have a constraint if A u <C 1, since 7 MK <C 1. 
This provides a lower bound on the number Nu of type II constraints; by choosing v from 
the m contracting directions, and summing over the remaining k < /i < v>, we obtain 

Nn > \m [n 2 - (m + 2) (n - \{m + 1))] . (5.16) 

But even if A u ^> 1 we can have a constraint, as long as A u '~f flK <ti 1. This depends on 
the particular problem at hand; hence, (5.16) is only a lower bound, but a fairly tight 
one for low dimensions. Table 2 enumerates the minimum number of type II constraints 
as a function of n and m. 

Note that when k, //, and v differ we can write 

1 

where the Lie bracket is 

d d 
, *A q : = E( 6 m)p q~ p (e v ), - E(^)p -q- p (e») q ■ (5-18) 

The type II constraints are thus forcing certain Lie brackets of the characteristic 
directions e a to vanish asymptotically. The geometrical implications of this, and perhaps 
a connexion to the Frobenius theorem [33] and the existence of submanifolds, remains 
to be explored. 
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Table 2. Lower bound on the number of type II constraints, as given by (5.16). The 
rows denote n, the columns the number of contracting directions m < n. 
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If we restrict to three dimensions, then there is at least one type II constraint; it 
can be written as 

©23i - @32i = ei • V x ei, (5.19) 

where Vo denotes a gradient with respect to the Lagrangian coordinates a. This 
constraint is the special case that was derived in Ref. [25]. 

5.3. Riemannian Curvature 

We now compare the approach of section 5.2 to the earlier attempts of Refs. [6,25], 
where the constraints were derived for two and three dimensional flows by examining 
the form of the Riemann curvature tensor associated with the metric g. We list the 
advantages of the present method. 

Nontrivial metrics can have curvature; a straightforward method of computing that 
tensor is through the use of Ricci rotation coefficients [34], 

d 

UJK^^^U^U^ — Ujy. (5.20) 



These satisfy the antisymmetry property io Kilv = —uj Kl/fl , and can be rewritten in terms 
of the $ of (4.8) as 

KIW - A' 1 . (5.21) 



v. 



In terms of the rotation coefficients, the Riemann curvature tensor is [34, p. 51] 

R^lVKa = J! ^ K ^oiiv — X! ^»<T W KHV 

i i 

— X! [^kt^ U aTU — U aT ^ LO KTV + U) KTa LO T ^ u — UJ UTK CO Tflu ] . (5.22) 

T 

If we use the relations (5.5) and (5.6) to solve for $ in terms of and ^, we can 
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rewrite (5.21) as 

Uk ^ v = OA a a ~~ *SW) + hl(Q Klll/ - Q^u) - A 2 K (Q flUK - 0^ K )| 

+^-6 liK V vlt -- j -6 VK V lu/ . (5.23) 

The form of the curvature obtained by inserting (5.23) into (5.22) is essentially the 
one obtained in three dimensions in Ref. [25]. In Ref. [6], the curvature was calculated 
directly from the Christoffel symbols. The constraints were then deduced by imposing 
the boundedness of the curvature tensor: some terms in (5.23) would appear to grow 
exponentially, so their coefficient must go to zero to maintain a finite curvature. In this 
manner, the type I and type II constraints were derived in two and three dimensions, 
backed by numerical evidence [6,25]. 

The approach used in the present paper to derive the constraints is advantageous 
in several ways: (i) It is valid in any number of dimensions; (ii) It avoids using the 
curvature, which is difficult to compute; (iii) There are no assumptions about the 
growth rate of individual terms in the curvature [25]; (iv) The convergence rates of 
the constraints are given explicitly [(5.10) and (5.15)]; (v) The number of constraints 
can be predicted [(5.12) and (5.16)]. The crux of the difference between the two 
approaches is that here we use the variable $ to estimate the asymptotic behaviour 
of the constraints directly, rather than relying on indirect evidence from the curvature. 
Thus, the derivation of the time-asymptotic form of $ is essential. 



6. Numerical Computations using the QR Method 

6.1. Basic Method 

The QR method, like the SVD method, avoids the numerical problems associated with 
evolving the Jacobian matrix M by using a judicious matrix decomposition. The QR 
decomposition says that any matrix, and in particular M, can be written 

M = QR, (6.1) 

where Q is an orthogonal matrix and R is upper-triangular. For our case, R has 
positive diagonal elements. The QR decomposition method of finding Lyapunov 
exponents is also called the continuous Gram-Schmidt orthonormalisation method by 
some authors [24,35], referring to the matrix Q being obtained from M by the Gram- 
Schmidt method. The QR method is an approximate version of the SVD method. The 
matrix Q is analogous to U in that it embodies the Eulerian information about the 
orientation of the ellipsoid (see section 4.1), and drops out of g as required. But the 
resulting expression g = R T R does not manifestly give a diagonalisation of g. Below 
in (6.7) and (6.8) we give the eigenvectors e CT and coefficients of expansion A CT in terms 
of R, though the expression is not exact but is exponentially accurate with time (see 
Appendix B). 



Derivatives and Constraints in Chaotic Flows 



16 



Let 

A^ '■— R^, Dfiv - = ^fx^ixvi = ' Dft^r^q. (6-2) 

That is, A is a vector containing the diagonal elements of R, D is a diagonal matrix 
with the A along the diagonal, and r is R with the /xth row rescaled by A M . The 
time-evolution of these quantities is [24, 26] 

= G^A,; (6.3) 

{-G^ // < iv; 

/i > !/; (6.4) 
\i = v\ 

9 A CT 

*imi= 12 ir- A ^ra q , H<q\ (6.5) 

cr=fl+l A 4 

where 

G := Q T GQ, A:=G + G T . (6.6) 

Equation (6.3) is identical to (4.2) for in the SVD method, and (6.4) is identical to 
the time-asymptotic form of (4.3) for U, given by (4.5). Hence, we expect A M and A M 
to have similar asymptotic behaviour, though their exact value differs. 

Unlike the SVD method, in the QR method the eigenvalues A^ and eigenvectors e M 
are not evolved directly. They can be recovered from the A M and the matrix r in the 
following manner. Let d be the lower-triangular matrix that effects the Gram-Schmidt 
orthonormalisation of r, that is 

W QH = 12 d »T r -rg » ( 6 - 7 ) 

T=l 

where = for ji < u, and W is orthogonal. The eigenvectors of g and corresponding 
coefficients of expansion are then 

(e,) q = W qit , A, = ^- , (6.8) 

to exponential accuracy with time (the relative error on e M and A 2 is of 
order (A M /A M _i) 2 ). Equation (6.8) is proved in Appendix B. 

By definition, the matrix W is obtained by Gram-Schmidt orthonormalisation 
of the upper-triangular matrix r. In performing this orthonormalisation, we have to 
compute the diagonal elements d^, so there is no extra work involved in correcting 
the A M if we are calculating the eigenvectors W. Note that this Gram-Schmidt procedure 
does not represent an extra overhead in solving the system of ODEs (6.3)-(6.5), as 
the orthonormalisation need only be effected at the end of the integration, when the 
eigenvectors are required. This orthonormalisation should not be confused with the 
continuous Gram-Schmidt orthonormalisation of the QR method, whose purpose is to 
evolve the orthogonal frame given by Q. 
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In most examples of applications of the QR method, the correction derived above 
to the eigenvalues is omitted [23,26]. The reason for this is that typically what is sought 
are the infinite-time Lyapunov exponents, 

A£°= Urn ^(logA^-logd^). 

Since the converge to constant values, they are irrelevant to the asymptotic value A^°. 
This means that it is possible to find the infinite-time Lyapunov exponents without 
solving for the eigenvectors. However, as mentioned in section 2, we are interested here 
in timescales much shorter than the convergence time of Lyapunov exponents, so we 
include the correction. 

We close this section by giving an explicit recurrence relation for the W qil and the 
lower-triangular matrix d: 



W, 



d 



fifj 



d 



E E r w w i° w p° 

(j=l q 



E^-ElE^wv 



(7=1 



-1/2 



= -d^ E E r itqW qa dau, n>v. 

u=v q 

These follow from the usual Gram-Schmidt procedure. 



(6.9) 
(6.10) 
(6.11) 



6.2. Lagrangian Derivatives 

We now proceed to obtain ordinary differential equations for the derivatives of the 
eigenvalues and eigenvectors of g, using the QR method. Define 



^KU ■= E W Q* 



<9 1ogA, 

dai ' 



4>k^u ■= E w<?« Q 



k/j, 



dai 



■■= E w qK 



dr 



hp 



q q 

where (f) Kllu = —<j) KV fi- The tensors ip and are the QR method analogues of \P and $ 
defined in (4.8) for the SVD method. The tensor £ has no analogue in the SVD method, 
but is used to obtain the Lagrangian derivatives of the eigenvectors W qix . 
Using the equations of motion (6.3)-(6.5) for A, Q, and r, we find 

i>™ = E(^ T W0- + E AnA- + Y KVU (6.12) 



<J 

^ ] A-CFfx (final* ^ ' A. av (pn^ia ^KVfl 1 

(7<fl (7>U 

L»q = Y,( WT WUUq 

(7 

q A CT r 



H<v, (6.13) 



<r=/x+l 



+ H<q. (6.14) 
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The driving term Y is 

n 

<T=K 

where 

■^VKfi ■ ^ ' Qku Qin Qtfi Qgi Qx^ 

is analogous to X of the SVD method, and is also symmetric in k and fx. The lower- 
triangular matrix d^ 1 = rW was defined in (6.7). 

In order to solve (6.12)-(6.14), we need to obtain the derivatives W T W. Because W 
is obtained from r via Gram-Schmidt orthonormalisation, the time derivatives of W are 
deduced from those of r by differentiation of (6.9). After multiplying that equation 
by W pu , with /i < u, we find 



(W T W), V = -d t 



E w pv r w + Y.( d ~ l )^ wT w\ 

P=fJ,+ l (7=1 



H<v. (6.15) 



Owing to the orthogonality of W, the matrix (W T W) flI/ is antisymmetric. 
Equation (6.15) defines a recurrence relation for the (W 1 ]^)^, starting with (W T W)i u , 
in terms of the time derivatives of r, given by (6.5). 

The recipe for finding the Lagrangian derivatives thus consists of solv- 
ing (2.1), (6.3)-(6.5) and (6.12)-(6.14) using a standard ODE integration scheme. In 
doing so, use must be made of the Gram-Schmidt procedure (6.9), which yields W and 
consequently also d via d~ x = r W. The matrix d must then be inserted into the re- 
currence relation (6.15) for W T W, allowing finally the full evaluation of the right-hand 
side of (6.12)-(6.14). The total number of ODEs involved is n(2n 2 + 3n + 3)/2; in two 
dimensions, this is 17, in three, 45. In evaluating the right-hand side of (6.12)-(6.14), 
the most expensive term to evaluate is Y, which scales as n 4 , obfuscating the cost of the 
Gram-Schmidt procedures for W and W T W . It is thus clear that this numerical method 
is not well suited to higher-dimensional dynamical systems. However, it is appropriate 
to applications such as chaotic mixing, where v is a two- or three-dimensional flow. 

We are not quite done yet: even though we can now solve the ODEs, they do not 
give the Lagrangian derivatives of the A M and W qil directly. The W qtl are obtained from 
the r vp via Gram-Schmidt orthonormalisation, so we need to proceed as we did for the 
time derivatives of W , (6.15), and take a Lagrangian derivative of (6.9). We obtain the 
recurrence relation 



n M— 1 

E wvCw + EtO^* 

p=/i+l cr=l 



H<v, (6.16) 



where 



p,q 



dW pv 
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is the analogue to in the SVD method. The recurrence relation is solved by first 
evaluating 9 K \ V and then incrementing //, always keeping \x < v. The antisymmetry 
of Kfll/ in fi and v means that we do not need to consider the // > v case. 

Finally, we need the Lagrangian derivative of in order to find the derivative 
of A^. Indeed, because of the correction to A M given in (6.8), we have 

where ^> KU was defined in (4.8), and 

d log d vv 



rinv ■= E W 



^]ku dp 



(6.17) 



<? 

is the correction. The explicit form for r\ is readily obtained in the same manner as (6.16) 
by differentiating (6.10), to yield 

n u—1 

P=U+1 <T=1 

Equation (6.17) is the same as (6.16) with /i — u, so that numerically both 9 and r\ can 
be obtained in the same loop. 

This completes the numerical procedure. As we mentioned in section 6.1, there 
is no real additional numerical burden involved in evaluating (6.16) and (6.17), as the 
Lagrangian derivatives of W and A M are not needed to solve the ODEs. These derivatives 
can be calculated as desired, either at regular intervals or at the end of the integration. 

There are two related but distinct numerical problems when finding the Lagrangian 
derivatives. The first is that the direction of fastest stretching of the flow dominates 
and must be isolated from the other directions, otherwise it quickly becomes impossible 
to extract subdominant directions because of lack of numerical precision. This is the 
problem we have solved with our method, by projecting along appropriate characteristic 
axes. The second numerical problem is that the exponentially growing quantities in the 
method eventually lead to numerical overflow (or underflow for exponentially decreasing 
quantities). In the QR method for the coefficients of expansion, (6.3) can easily be 
rewritten as an equation for log A M , replacing the exponential behaviour by linear growth 
(or decay) in time. The same cannot be done in (6.12)-(6.14) because the rescaling of 
the driving term introduces a large damping term that makes the system extremely 
stiff (because the rescaling itself is time-dependent). But overflowing only becomes a 
problem if we solve the system for very long times (on the timescales necessary for the 
Lyapunov exponents to converge). 

A final note on the stability of the algorithm: in Refs. [23] and [35] it is shown 
that the numerical integration of the orthogonal matrix Q is unstable (the matrix 
loses orthogonality), unless the full spectrum of eigenvalues is calculated, in which case 
the algorithm is neutrally stable. Since we always assume we are computing the full 
spectrum, the stability of Q is not a concern. 
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Figure 1. The type I constraint given by (5.13), for the ABC flow with A = B = 5, 

C = 2, computed with the direct method of section 3 ( ) and with the QR method 

of section 6 (- ) . The direct method becomes unreliable after t ~ 7 due to roundoff 
error. The QR method unambiguously exhibits the convergence of the constraint to 
zero, which agress well with the predicted convergence rate A 3 ( ). 

6. 3. Numerical Verification of Constraints 

In figure 1, the type I constraint given by (5.13) is shown for the ABC flow with A = 
B — 5, C — 2 [11]. (These parameter values give a large chaotic region and make the 
convergence of the constraints faster, but the constraints are also satisfied for the more 
usual values A — B — C — 1.) The constraint is computed with the direct method of 
section 3 and with the QR method of section 6. It would be difficult to make a case 
for the constraint converging to zero based on the direct method: the noise starting 
at t ~ 7 reflects the effects of limited numerical precision inherent to the method as the 
elements of M become exponentially large. The QR method, however, has the constraint 
reaching 10 -12 before precision problems set in (this is not a flaw in the method: the 
terms in (5.13) cannot cancel beyond the number of digits of precision represented by 
the machine). The constraint is predicted by (5.13) to converge as A 3 , which is also 
shown in figure 1. 

Figure 2 shows a plot of the type II constraint given by equation (5.19), also for 
the ABC flow. The constraint converges as (A 2 /Ai) 2 , as predicted by (5.15). The same 
comments as for figure 1 apply regarding numerical precision. 




Figure 2. The type II constraint |ei • Vo x ei| = I6231 — 6321! [equation (5.19)] for 
the ABC flow with A = B = 5, C = 2, computed with the direct method of section 3 

( ) and with the QR method of section 6 ( ). The direct method becomes 

unreliable after t ~ 7 due to roundoff error. The QR method clearly illustrates the 
convergence of |ei • Vo x ei| to zero, at the predicted rate of (A 2 /Ai) 2 ( ). 



7. Discussion 

Lagrangian coordinates can greatly simplify the form of partial differential equations, 
specifically equations of an advective nature. We have taken the viewpoint that to fully 
characterise quantities expressed in the Lagrangian frame it is necessary to know how 
to compute derivatives with respect to these Lagrangian coordinates. This amounts 
to understanding how the Lagrangian frame itself (as defined by the characteristic 
eigenvectors and the coefficients of expansion) vary under small changes of initial 
conditions. Obtaining such derivatives with accuracy is difficult in chaotic flows because 
the stretching rates of phase-fluid elements vary greatly along different directions. 

The Lagrangian derivatives can be computed by differentiating existing methods for 
finding Lyapunov exponents and eigenvectors. Direct differentiation of the equations of 
motion is useful only for short times. For long times limited numerical precision becomes 
problematic and a decomposition method is needed. The SVD method proved useful in 
deriving the asymptotic form of the Lagrangian derivatives and in deriving differential 
constraints. The singularities possessed by the SVD method and its large number of 
components make its use in numerical computations difficult. The QR decomposition 
method is more appropriate for numerical implementations, but is less transparent 
than the SVD method. We used the QR method to accurately verify the differential 
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constraints derived in the paper. 

The techniques described here apply only to the first derivatives of the various 
quantities, which actually depend on the second derivatives of the vector field. First 
derivatives are sufficient for the study of many systems, including the advection- 
diffusion equation and the dynamo problem. This seems paradoxical because both 
the advection-diffusion equation and the resistive magnetic induction equation involve 
second derivatives in space, but when transformed to Lagrangian coordinates the 
equation only involves first derivatives of the metric tensor g [6, 12]: the second 
derivatives are only applied to the initial data. 

The results derived in the Lagrangian frame can easily be adapted to the Eulerian 
picture with only minor modifications. The crux of the difference lies in holding the 
Eulerian final condition fixed, and integrating the initial condition backwards in time. 
The translation between the Eulerian and Lagrangian pictures is outlined in Appendix 
C. 

Some of the differential constraints derived in this and earlier papers have been 
applied to the study of the advection-diffusion equation [6, 12, 31]. In particular, 
in Ref. [12] a type I constraint (Section 5.2.1) is used to obtain an effective one- 
dimensional diffusion equation for chaotic flows. In Ref. [36] the Eulerian form of the 
type I constraint, equation (C.5), is used to derive a power-law relationship between 
the curvature of a material line and the amount of stretching the line has undergone. 
This same constraint is used in [8,9] to derive an invariant measure of the spatial length 
distribution of material lines in two dimensions. Finally, in Ref. [37] a type II constraint 
(Section 5.2.2) is used to show that the onset of dissipation in the kinematic dynamo 
occurs much later than a straightforward estimate indicates. This is because the leading- 
order behaviour of the power dissipation (Ohmic heating) in the dynamo is proportional 
to a type II constraint, so it does not grow as fast as expected. 

The behaviour of second and higher derivatives has not been investigated. Whilst 
in principle the method could be extended to cover such cases, the complexity of 
the calculation and the smoothness requirements on v are prohibitive. A study of 
the consequences of the degeneracy of Lyapunov exponents — as occurs for instance in 
Hamiltonian systems — remains to be done. 
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Appendix A. Asymptotic Behaviour of the Lagrangian Derivatives 

In this appendix we use the equations of motion (4.9)-(4.11) to derive the asymptotic 
behaviour (4.13)-(4.15) of the Lagrangian derivatives, as defined by (4.8). 
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For /j, < v, assuming A M 3> A u , the last term in (4.10) is 
F) Fir* fir 1 

E v 9K £(u T uu = - E v qK ^ - ^ ^ + 2^ v - *«0- (A.i) 

The first term in (A.I) is 

E K« = a k + e ^ + E W (A.2) 

q a a 

This shows that the second term in (A.I) can be neglected compared to the first because 
it is smaller by a factor ^ 2 . We also neglect the third term, which couples $ and \l/ (it 
is straightforward to go back and check that the neglect is justified). 
After these approximations, the evolution equation for & KfMU is 

a 

+ 52(V T V) ffK <f> aiJH , - A K X VKIi , (A.3) 

a 

for fj, < v. This is a linear system in $, with nonconstant coefficients and a driving term 
given by — A K X L , K ^ l . The only term that couples $ K ^'s with differing k is the V T V one, 
which is small compared to U T U. We neglect this term (again, after the fact it is easy 
to check that the neglect is consistent), and rearrange (A.3) to give 

$K»V = (Gvv ~ G^KIXV + E A va $ KtuJ - E ~ A « ( A - 4 ) 

cr>u cr<fl 

Let us ignore the driving term for now and consider the homogeneous solution ^ fil/ . 
Through a judicious ordering of the § Kflu that gives the linear part of (A. 4) a triangular 
structure, it can be shown that ~ 7 MJ y- If we assume that the motion of the system 
takes place in a bounded region of phase space, X UKfl is also bounded (we also assume 
that v is twice differentiable and that its second derivative is Lipschitz). Then the 
inhomogeneous driving term k K X UKfl asymptotically goes as A K . (A similar argument 
was used in section 4.1 to show convergence of V.) Asymptotically, then, in (A. 4) either 
the exponentially decaying linear part or the driving term dominates, depending on 
which has the larger growth rate. We conclude that 

$ KIU/ = max (A K , 7^) t > 1, (A.5) 

where § Kflu is some function that neither grows nor decays exponentially. 

Next, we investigate the asymptotic behaviour of Its time evolution is given 
by (4.9), which after inserting our asymptotic solution for $ becomes 

= £(^U + max (A K , lvv ) 

where \I/^! ve is some non-exponential function. Notice that ^^'s with differing v are 
uncoupled. The matrix V T V has elements that are decreasing exponentially, so we can 
solve the system perturbatively. The solution, valid to first-order in V T V, is of the form 

^ = max(A K ,7 KK , 7w )^ + ^~, t > 1. (A.6) 
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where the time dependence of ty KV is non-exponential. 

Finally, having derived an asymptotic form for $ and we can do the same 
for the Lagrangian derivatives of V, as embodied by [Equation (4.8)]. For /j, < u, 
assuming A M ^> A„, we write (4.11) for in the approximate form 

<sw = E [(v T v)„ K + (v T v)^ e Kau - (v T v) au e Ka , 

+^ I1U max (A (A.7) 

where 0^ c is a term with possible time dependence but without exponential behaviour. 
The matrix V T V, given in (4.5), becomes exponentially small with time. We can thus 
solve (A.7) perturbatively, yielding 

e KfJ/l/ = max (7 M1 ,A K , 7 KK , 7^, ~j vv ) K/li , + 0~„, t > 1, (A.8) 

where the time dependence of Kflu is non-exponential. 

A comment about the perturbation expansions in small V T V used to obtain (A. 6) 
and (A.8) is in order. For a small parameter e, a perturbative expansion solution to an 
equation of the form 

y = e ay + (3 

is valid only for et <^ 1. So for large time, at fixed e, the solution must eventually 
become invalid. However, in our case the parameter e, corresponding to V T V, actually 
decreases exponentially in time (for nondegenerate eigenvalues A M ). Thus, et ^1 even 
for large t, and the expansion remains valid. 



Appendix B. The Eigenvalues of g and the QR Method 

In this appendix we prove that (6.8) gives the asymptotically correct value of the 
eigenvectors e M and coefficients of expansion A^ (the square root of the eigenvalues 
of g) . This result was shown in Ref . [24] . We present here a different proof, proceeding 
by induction and deriving the eigenvectors and eigenvalues together. 

Let d be the lower-triangular matrix that performs the Gram-Schmidt 
orthonormalisation of r, that is 

W T = dr, (B.l) 

where d is lower-triangular and W is orthogonal (this is simply (6.7) in matrix form). 
The matrix d is nonsingular, so we can invert (6.7) and write r = d~ 1 W T . 

Using the definition of the metric, g = M T M, and the QR decomposition (6.1), we 

have 

g W = W3 r , where 5" := (d^f D 2 (cT 1 ). (B.2) 

If jF were diagonal, then we would be done; the matrix 5F is not diagonal, but we show 
that it can be made so by exponentially small corrections of the W. 
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Writing out 3 explicitly, we find 

n 

^ = E (O^X (O™- 

cr=ma,x(fi,u) 

The max(//, v) is the lower bound of the sum owing to the lower-triangular form 
of d. Looking at (6.5), it is clear the the r converge to constant values, because 
for o > fi, A a /A tl — > exponentially in time. Hence, by their definition, W and d also 
converge to constant values in time. All of the exponential behaviour is thus embodied 
in A CT . Keeping only the dominant term, we have 

Note that since d is triangular we have = d~*. For the first column of W , it 

follows that 

A 2 

E 9pq Wql = Wl ~ ^ W Pl > 
g cr "11 

showing that the first column of W is an eigenvector of g with eigenvalue (Ax/rfn) 2 . 
Let 

= W v» - E 7^ ^ ^ (B.3) 

which represents an exponentially small correction to the matrix element W PfM , 
because 3^ ~ A 2 A 2 for a < /i. Assume that the columns W' qv , v < /i, are 
eigenvectors of g with eigenvalue (A u /d uu ) 2 . We show by induction that with the (small) 
correction (B.3) the column W is an eigenvector of g with eigenvalue (A^/rf^) 2 . 
Using 5F = W T gW and (B.3), the corrected matrix element 3^, with v < (j,, is 

^ = ^ a* wv = v - E % ^ ^ov, «/ < /i- 

(7=1 ^cr 

We use the induction hypothesis that J av = (A v /d vv ) 2 8 au when both v and a are less 
than fi, and find 

fi-l rf 2 A 2 

to exponential accuracy. The corrected diagonal element 5F' is 

M- 1 Al M-l M-l J2 J2 

qr/ _ qr _ o \ ^ era qr qr _i_ \ ^ \ ^ "'erg- "tt qr qr qr 

J MM ~~ A 2 ^ + Z^Zj A2 A2 ^ T °" r 

cr=l ^cr cr=l r=l ^cr 

M-l j2 
(T=l ^cr 

to leading order. Thus, to exponential accuracy W 7 ^ is indeed an eigenvector of g with 
eigenvalue 3^ = (A M /c^) 2 . 
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To complete the proof, we need to show that the columns W qu , with v > /j,, are not 
modified by the correction (B.3). We have 

d 2 

T' = T - V — T T v > u 

which to leading order is 

%u = (d-'U {d- v U K + o (A- A 2 ) ~ ^, 

showing that the correction can be neglected. 

We have thus proved by induction that the correction (B.3) to W makes 5F diagonal 
to leading order, leaving its diagonal elements unaffected. However, the correction (B.3) 
is of order (A^/A^i) 2 , which is exponentially small with time. We conclude that the 
eigenvectors of g and corresponding coefficients of expansion are 

(e„), = W q „, A M = ^ , 

to exponential accuracy with time. The relative error on e M and A 2 is of 
order (A^/A M _i) 2 , as can be seen from the leading-order correction in (B.3) and (B.4). 



Appendix C. The Eulerian Perspective 

The results derived in the paper regarding Lagrangian derivatives can readily be adapted 
to an Eulerian framework. The Lagrangian derivatives can be regarded as measuring 
the effect of an infinitesimal change in the initial condition of a trajectory. Conversely, 
one can regard Eulerian derivatives as the effect of an infinitesimal change in the 
final condition of a trajectory, with the integration being performed backwards in 
time. This is the viewpoint taken in studies of the alignment of material lines with 
the Eulerian unstable manifold of a system, a phenomenon referred to as asymptotic 
directionality [8,9,22,38,39]. 

In our framework, the Eulerian characteristic directions are computed from the 
metric 

n 

hV(t,t ,x):=J2M\M\, (C.l) 
p=i 

or h := MM T = g T . We have explicitly written the dependence on initial time 
because we are now interested in evolving £ towards — oo, whilst holding the Eulerian 
coordinates t and x fixed. The dynamical system (2.1) and the SVD equations (4.2)- 
(4.4) are also evolved backwards in time: the method yields M _1 , so that we must take 
the inverse of the resulting A CT to obtain the "forward-time" coefficients of expansion. We 
assume the forward-time coefficients have then been reordered in the usual decreasing 
manner, so that A 1 is still the fastest-growing coefficient (the columns of U and V 
are also reordered). The asymptotic behaviour of the Eulerian derivatives will thus be 
the inverse of their Lagrangian counterpart. The columns of V now contain vectors 



Derivatives and Constraints in Chaotic Flows 



27 



associated with the Eulerian frame. The relevant Eulerian definitions corresponding 
to (4.8) simply involve replacing d/da by d/dx to reflect the fact the the derivatives 
are now taken with respect to the Eulerian coordinates. Their asymptotic behaviour is 

e^ = max(A^,7^)e^, (C.2) 

* = max (A" 1 , 7kk , 7w ) *l + , (C.3) 

= max (t^A" 1 , 7kk , 7 W , iJj ^ + (C.4) 

where the £ superscript reminds us that these are Eulerian quantities. The discussion of 
the Hessian and constraints in section 5 is the same for the Eulerian derivatives, except 
that the coefficients of expansion corresponding to contracting directions are replaced 
by the inverse of those of the expanding directions. For instance, the type I constraint 
given by (5.13) becomes 

E \h\ 1/2 ^ 0r V2 (ef )i) + E(ef )^ log Ai - max (Ar 1 , 7l 2 x ) - 0. (C.5) 

i i 

This is a more general form of the relation derived for the 2D incompressible case in 
Refs. [8,9], used to derive an invariant measure of the spatial length distribution of 
material lines. 
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